Eccentric-disk models for the nucleus of M31 

Hiranya V. Peiris^ and Scott Tremaine^ 
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 



(N 



^- ABSTRACT 

O 
O 

CN ■ We construct dynamical models of the "double" nucleus of M31 in which the 

nucleus consists of an eccentric disk of stars orbiting a central black hole. The 
principal approximation in these models is that the disk stars travel in a Kepler 
potential, i.e., we neglect the mass of the disk relative to the black hole. We 
consider both "aligned" models, in which the eccentric disk lies in the plane of 
^ . the large-scale M31 disk, and "non-aligned" models, in which the orientation 

^ . of the eccentric disk is fitted to the data. Both types of model can reproduce 

. the double structure and overall morphology seen in Hubble Space Telescope 

^ . photometry. In comparison with the best available ground-based spectroscopy, 

. the models reproduce the asymmetric rotation curve, the peak height of the 

dispersion profile, and the qualitative behavior of the Gauss-Hermite coefficients 
'q_|! hs and h^. Aligned models fail to reproduce the observation that the surface 

brightness at PI is higher than at P2 and yield significantly poorer fits to the 
kinematics; thus we favor non-aligned models. Eccentric-disk models fitted to 



O 
H 



c§ ! ground-based spectroscopy are used to predict the kinematics observed at much 

higher resolution by the STIS instrument on the Hubble Space Telescope (Bender 
^ ! et al. 2003), and we find generally satisfactory agreement. 

Subject headings: galaxies: individual (M31) — galaxies: nuclei — galaxies: 
photometry — galaxies: kinematics and dynamics 



1. Introduction 

The curious asymmetric nucleus of the Andromeda Galaxy (M31=NGC 224) has in- 
trigued astronomers since it was first resolved by the balloon-borne telescope Stratoscope 
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II in 1971 (Light et al. 1974). Thus the nucleus of M31 was a prime target for the Hubble 
Space Telescope (HST). Early HST photometry (Lauer et al. 1993, hereafter L93; King, 
Stanford & Crane 1995) revealed that the apparent asymmetry arose because the nucleus 
was double, with two components separated by O'.'S; the fainter peak P2 corresponds closely 
to the dynamical and photometric center of the bulge, and the brighter peak PI is off-center, 
with the nucleus confined roughly within 2" of P2. Ground-based spectroscopy, although 
unable to resolve the two components, revealed a prominent velocity-dispersion peak and 
strong rotation-curve gradients that implied the presence of a massive black hole (Dressier 
& Richstone 1988; Kormendy 1988; Bacon et al. 1994; van der Marel et al. 1994). 

It was natural to assume at first that PI and P2 were orbiting star clusters. However, 
the orbit of these clusters would decay by dynamical friction from the bulge in < 10^ yr, so 
the present configuration is highly improbable. Difficulties with this and other hypotheses 
are discussed by L93 and Tremaine (1995, hereafter T95). 

A more attractive possibility, suggested by T95, is that the nucleus consists of an ec- 
centric disk of stars orbiting a massive black hole (hereafter BH). In this model, the stars 
travel on eccentric orbits with approximately aligned apsides; the BH is located at P2, and 
the bright off-center source PI is the apoapsis region of the eccentric disk. T95 argued 
that an eccentric disk can reproduce most of the features seen in the HST photometry and 
ground-based spectroscopy. 

The model in T95 had several serious limitations. For example, it uses only three Keple- 
rian ringlets rather than a continuous distribution of orbits; the eccentricity and inclination 

dispersion in the disk are modeled approximately using a Gaussian point-spread function 
(PSF); the effect of the self-gravity on the disk on the stellar orbits is neglected; and the 
model parameters are chosen by eye, rather than by numerical parameter fitting. 

The eccentric-disk model requires that the apsides of the disk stars precess uniformly 
so that the disk maintains its apsidal alignment. The uniform precession rate is simply 
the pattern speed of the eccentric disk. Statler (1999) investigated the restrictions that 
this condition imposes on the parameters of the eccentric disk. He argued that uniform 
precession requires that the disk have a steep eccentricity gradient and a change in direction 
of the eccentricity vector, in the sense that stars in the inner part of the disk have apoapsides 
aligned on the PI side, while stars in the outer part have periapsides aligned on the PI side. 

Salow & Statler (2001) and Sambhus & Sridhar (2002) have presented eccentric- disk 
models for the nucleus of M31 in which the apsides of the disk stars precess uniformly. 
These models are constructed by numerical integration of orbits in the combined potential 
of a BH and a nuclear disk that is determined self-consistently from the mass distribution of 
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the orbits; in this respect their models are better than the ones in this paper, which neglect 
the effect of the disk potential on the orbits. However, both papers assumed that the disk is 
razor-thin, which is unlikely to be correct — two-body relaxation alone will thicken the disk 
to an axis ratio ~ 0.2 if the disk age is comparable to the Hubble time (see §4.2). A^-body 
simulations of the disk by Bacon et al. (2001, hereafter BOl) were also mostly restricted 
to two dimensions; some three-dimensional simulations were carried out but not described 
in detail. The assumption of a razor-thin disk drove Sambhus & Sridhar (2002) to models 
in which the eccentric disk is not aligned with the larger M31 disk (we shall confirm the 
claim of T95 that a thick eccentric disk ahgned with the M31 disk can provide a reasonably 
good fit to the data, but hke Sambus & Sridhar and BOl we find that non-aligned models 
fit better than aligned ones). Salow & Statler (2001) did not have to consider a non-aligned 
disk because they compared the photometry to the data only along the major axis. Salow 
& Statler estimate a pattern speed of 14 km s^^ pc~^ while Sambhus & Sridhar (2002) find 
16 km s^^ pc^^. The close agreement between these two estimates does not imply that they 
are reliable, since the pattern speed makes only a small contribution to the overall kinematics 
(at the P1-P2 separation a pattern speed of 15 km s~^ pc~^ corresponds to a velocity < 20% 
of the peak rotation velocity), and since the dynamical models used in the two papers are 
quite different (the mass ratio of the eccentric disk to the BH is 0.03 in Salow & Statler and 
0.65 in Sambhus & Sridhar). 

BOl carry out A^-body simulations of an eccentric disk that is intended to resemble 
the nucleus of M31. Their experiments demonstrate that long-hved eccentric configurations 
arise naturally from a variety of non-axisymmetric (lopsided) initial conditions. Their best- 
fit simulation reproduces the photometry and rotation curve of the M31 nucleus reasonably 
well. However, once again they assume that the disk is thin (axis ratio 0.1). They estimate 
that the pattern speed is ~ 3 km s~^ pc~^, far smaller than the estimates of Salow & Statler 
and Sambhus & Sridhar. 

A variety of new, high- resolution, space- and ground-based photometry and spectroscopy 
of the nucleus of M31 has become available since 1995. These include multicolor WFPC-2 
images from HST with greater signal-to-noise ratio (S/N), spatial resolution, and dynamic 
range than L93 (Lauer et al. 1998, hereafter L98), which have been analyzed in detail by 
Peng (2002); near- infrared (JHK) photometry using both adaptive optics (Davidge et al. 
1997) and the NICMOS instrument on HST (Corbin, O'Neil & Rieke 2001); ground-based 
long-sht spectroscopy with slit width=0'.'35 and seeing FWHM=0'.'64 (Kormendy & Bender 
1999, hereafter KB99); spectroscopy with the Faint Object Camera (Statler et al. 1999, 
hereafter S99) and Faint Object Spectrograph of HST (Tsvetanov et al. 1998); integral-field 
spectroscopy using adaptive optics, with FWHM=0'.'5 (BOl); and long-slit spectroscopy by 
the Space Telescope Imaging Spectrograph (STIS) on HST (BOl, Bender et al. 2003). 
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Given the wealth of new observations and the hmitations of existing models, the time 
is ripe to revisit the eccentric-disk hypothesis and confront more realistic models with more 
accurate, varied and higher resolution data. In §1.1 we review the rich phenomenology of the 
M31 nucleus, and in §1.2 we briefly describe the eccentric-disk model. In §2, we construct 
an improved kinematic eccentric-disk model with a range of free parameters that wc can fit 
to the observations. In §3 we describe the details of the numerical implementation of the 
model and parameter fitting. We describe our results in §4, and provide a discussion in §5. 

We shall not attempt to fit or discuss all of the observations described above. In- 
stead, wc shall focus on fitting the HST photometry described by L98 and the ground-based 
spectroscopy of KB99. As a test of our model, we shall compare its predictions to HST 
spectroscopy from Bender et al. (2003) (see BOX for an earher analysis of the same data). 

Constructing general, self-consistent models of stellar systems to compare with obser- 
vations is a difficult task that is important in many areas of galactic structure. The state 
of the art is orbit-based models of axisymmetric systems, in which the distribution function 
can depend on up to three integrals of motion (Cretton et al. 1999; Gebhardt et al. 2000; 
Verolme et al. 2002). Eccentric-disk models are substantially more comphcated, because 
they are non- axisymmetric and involve up to five integrals of motion. To keep the problem 
manageable we shall adopt a number of simplifying approximations. Of these, the most 
important is that the stars follow Kepler orbits in the gravitational field of the BH; that is, 
we neglect the effect of the self-gravity of the eccentric disk on the stellar orbits. The effects 
of this approximation are discussed briefly in §5. 

Throughout this paper, we adopt KB99's estimates of the distance and foreground 
extinction to M31, d ^ 0.77 Mpc and Ay = 0.24 (Burstein & Heiles 1984). At this distance 
1" = 3.73 pc. 

1.1. Observations 

The observational facts can be summarized as follows: 

1. The nucleus is composed of two components, separated by 0'.'49 and labeled PI and P2. 

The two components cannot be decomposed into a superposition of two systems with 
elliptical isophotes (L93). The total bulge-subtracted magnitude of both components 
of the nucleus isV = 12.55 ± 0.2, corresponding to a total luminosity L = 6.0 x IO^Lq 
(KB99). The same double structure is seen in near-infrared images, demonstrating 
that the apparent duplicity of the nucleus is not due to a dust band (Davidge et al. 
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1997; Corbin, O'Neil & Rieke 2001). There is no evidence that the photometry is 
affected by spatially varying obscuration. 

2. The position angle of the hne joining PI and P2 is 43°±1° (as usual, measured eastward 
from north), while the position angle of the isophotes within I'/O of the center of PI is 
about 63° (L93); thus there is an isophote twist of about 20°. The position angle of the 
outer parts of the nucleus (1-2") is 55° ±1°, and the position angle of the large-scale 
M31 disk is 38°. 

3. A compact blue source is centered on P2 (we shall call this P2B); as a result of this 
source, P2 is actually brighter than PI in both the near- and far-UV (L98; King et al. 
1995). King et al. (1995) suggest that P2B is a low-level AGN; however, the higher S/N 
and resolution of L98's data reveal that P2B is extended, with a half-power radius 0'.'2 
and position angle 62° ± 8° (BOl). STIS spectra reveal that P2B is a cluster of carly- 
typc stars, with the remarkably high velocity dispersion of 800-950 km s~^ (Kormendy, 
Bender & Bower 2002; Bender et al. 2003). 

4. The stellar component of P2, which surrounds P2B, exhibits a weak cusp, I{r) oc r'"^ 
with 7 ~ 0.1 (L93). The peak in the stellar component appears to be slightly offset 
from P2B; the measured offset is strongly wavelength-dependent because of the color 
difference. 

5. The central surface brightness of PI is /xy = 13.4 mag arcsec"^ (as measured in a slit 
0'/22 wide; see L93) and the major-axis core radius is about 0'/4. In contrast to P2, 
PI exhibits no compact blue source and no cusp in the starlight (L93; King et al. 
1995); the central surface brightness of P2 in the same slit is 13.7 mag arcsec"^, 0.3 
mag fainter. 

6. The source P2 is very close to the center of the galaxy (within O'/l according to L98) as 
defined by the centroids of the bulge isophotes just outside the nucleus; KB99 estimate 
that P2B is displaced from the bulge center by 0'.'07 in the anti-Pl direction. 

7. The V — I color of the nucleus appears to differ from the surrounding bulge (although 
L98 and BOl disagree on the color difference); the color difference implies that the 
stellar populations in the bulge and nucleus are different, while the populations in 
PI and P2 (outside P2B) are the same (L98). KB99 reach a similar conclusion from 
absorption-line strengths. 

8. The gradients in both the rotation curve and the velocity-dispersion profile are unre- 
solved or marginally resolved, that is, they become steeper as the resolution improves. 
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Both profiles are asymmetric. KB99 find that the maximum (bulge-subtracted) rota- 
tion speed is — 236±4km s~^ on the anti-Pl side but only 179±2 km s~^ on the PI side 
(see Tabic 3). S99 claim that they have resolved the central gradient in the rotation 
curve, at 1000km s~^ arcsec^^; the rotation curves from S99 and KB99 are consistent 
when S99's HST results are smoothed to KB99's ground-based resolution. The peak 
velocity dispersion is displaced from P2 in the anti-Pl direction, by O'.'OG ± O'.'OS ac- 
cording to S99 (without bulge subtraction) or 0'/13 according to KB99 (with bulge 
subtraction). The peak dispersion is 287 ± 9km s~^ according to KB99 (with bulge 
subtraction), and 440 ±70 km s~^ according to S99 (without bulge subtraction). With 
bulge subtraction the peak dispersion in the S99 data would be even higher, approach- 
ing a factor of two larger than KB99. Some of this difference presumably refiects the 
higher resolution of HST but it is also possible that the S99 dispersions are system- 
atically high. Remarkably, at 1" away from the nucleus (in either direction along the 
major axis) the dispersion has fallen by a factor of three or more, to < 100km s^^. 
BOl argue that S99 have made a 0'.'074 error in registering the KB99 spectroscopic 
data with their HST data. 

9. The zero-point of the rotation curve (relative to the systemic velocity of the bulge) is 
displaced from P2B toward PI. KB99 find that the displacement is '.'051 ±0 '.'014 S99 
find '.'13±0 '.'05 (composed of '.'16±0 '.'05 offset from the steUar peak of P2 minus '.'025 
offset between the stellar peak of P2 and P2B) and BOl find 0'/031. These differences 
may reflect the small offset between P2B and the stellar peak of P2, which means that 
the measured position of the peak brightness of P2B is wavelength-dependent, as well 
as differences in spatial resolution (see BOl for discussion). 

10. KB99 also measure the Gauss-Hermite coefficients of the line-of-sight velocity distribu- 
tion. They find that the zero-point of hs is displaced by about 0'/04 from the zero-point 
of the rotation curve, in the anti-Pl direction. 

1.2. The eccentric-disk model 

In this model, the nucleus contains a single BH located at the center of P2B. This is 
consistent with (i) the presence of a compact cluster of young stars, perhaps formed from 
a gas disk surrounding the BH (the stellar collision rate is too slow to explain the presence 
of these stars; see Yu 2003); (ii) the high velocity dispersion of the P2B cluster; (ni) the 
weak surface-brightness cusp observed in P2 outside P2B; (iv) the presence of an unresolved 
dispersion peak in the old stars close to P2 but outside P2B. 
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Most or all of the stars in the nucleus are assumed to lie in a disk surrounding the BH. 
The mass of the stellar disk is assumed to be small compared to the mass of the BH (see §5 
for a discussion of this approximation). Thus, the stars travel on approximately Keplerian 
orbits, and the asymmetry of the nucleus arises because the orbits are eccentric and the 
mean eccentricity vector is non-zero (i.e. the apsides are approximately aligned). The off- 
center source PI marks the apoapsis region of the disk, which is bright because stars linger 
at apoapsis, while P2 (outside P2B) represents a combination of disk stars with smaller 
semimajor axes and stars with periapsis near P2 and apoapsis near PI. 

The disk is assumed to be in a steady state (there is presumably slow figure rotation as 
the apsides precess due to the self-gravity of the disk and the tidal field from the bulge, but 
we neglect the effect of figure rotation on the kinematics) . The center of mass of the stellar 
disk plus the BH must therefore lie at the center of the bulge. It is natural to assume that 
the age of the disk is comparable to the age of the galaxy, i.e. ~ 10^*^ yr, although this is 
not required for the model. 

The eccentric-disk model is consistent with the color and fine-strength observations, 
which imply that the stellar populations in the bulge and nucleus are different, while the 
populations in PI and P2 (outside P2B) are the same. It is also consistent with the obser- 
vation that the BH at P2 is located close to the center of the bulge, which is expected if the 
BH mass is much larger than the mass of the stars in the eccentric disk. Alternative models 
in which both both PI and P2 are stellar clusters bound to BHs would suggest that both 
should have a cusp and perhaps a compact blue source, but PI has neither. 

It is natural to assume at first that the eccentric nuclear disk lies in the plane of the 
large-scale M31 disk ("aligned models"), although we shall find that "non-aligned" models 
in which this is not the case provide better fits to the data. In aligned models the disk must 
be relatively thick: M31 is highly inclined (inclination 77°), so the disk is seen nearly edge-on 
and the isophotcs of a thin disk would be fiattened too strongly. 

The T95 model, although not unique, correctly predicted several features that were 
discovered in subsequent observations, such as the displacement of the zero-point of the 
rotation curve from P2B toward PI, and the shape and amplitude of the asymmetry in the 
rotation curve. S99 explored a number of similar eccentric-disk models and were able to 
improve the kinematic fit, but found it difficult to reproduce simultaneously the location 
of the zero-point of the rotation curve, the steepness of the rotation curve, and the overall 
photometric symmetry of the nucleus about P2 beyond the distance of PI (see also the 
discussion of S99's results in §§2.3 and 5.3 of BOl). 

In order to gain some intuition for the features of a eccentric-disk model, we review here 
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the simple case of an infinitesimally thin, cold disk with a BH at the origin (Statler 1999). 
Stars are assumed to orbit on confocal, nested ellipses. It is straightforward to show, using 
standard formulae for Kepler orbits, that the surface density along the x-axis is given by 

E(a)^^^ , (1) 

^ ^ 27ra ^/r^(l ± e ± ae') ^ ' 

where the + or — sign is chosen at apoapsis or periapsis respectively, \i{a) = dm/ da, dm is 
the mass contained between a and a + da, and e' = de/da. 

Denoting the surface brightness at periapsis and apoapsis as Ep(a) and Sa(a) respec- 
tively, the brightness ratio at a given value of a is 

S„(a) _ l-aeV(l-e) 



Ep(a) l + ae'/{l + e)' 

To meet the condition that the disk have a substantially higher surface brightness at apoapsis 
we require e' < and either e or \ae'\ of order unity. Therefore, as already noted by T95 
and S99, successful models require a strong negative eccentricity gradient in the interval of 
semimajor axes corresponding to PI (see Fig. 1). 



2. A model for an eccentric disk 

2.1. Coordinate systems 

We start by establishing our notation and defining the coordinate systems that we will 
employ. We use the BH as the origin of all of our coordinate systems, and assume that the 
BH coincides with the center of the blue source P2B. 

Our first coordinate system, the "sky-plane" system, is denoted by {X,Y,Z). The 
(X, Y) plane is the sky plane; the positive X-axis points west, the positive F-axis points 
north, and the positive Z-axis points along the line of sight toward the observer. In this 
system the line-of-sight velocity is Vios = —Z; the minus sign is necessary for consistency 
with the usual convention that objects receding from the observer have positive velocity. Our 
second coordinate system is the "disk-plane" system, denoted by {x, y, z). The {x, y) plane is 
defined to be the symmetry plane of the eccentric nuclear disk. The x-axis is aligned with the 
major axis of the nuclear disk (more precisely, the distribution of eccentricity vectors [defined 
at the end of §2.2] of the disk stars is assumed to be symmetric about the x-axis). A third 
"orbital-plane" coordinate system {x',y',z'), different for each star, is defined so that the 
{x', y') plane is the orbital plane of the star, the positive x'-axis points toward the periapsis 
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of the star, and the positive z'-axis is parallel to the star's orbital angular-momentum vector. 
All three coordinate systems are right handed. 

The transformation between the orbital-plane and disk-plane coordinates is given by: 

1 
COS / — sin / 
sin I cos / 

' ■ '(3) 

where / is the inchnation of the orbit relative to the disk plane, Q is the longitude of the 
ascending node of the orbit on the disk plane (the angle in the disk equatorial plane from 
the X-axis to the ascending node, where the orbit has z — 0, z > 0), and u is the argument 
of periapsis (the angle in the orbital plane from the ascending node to the x' axis). 

The transformation between the disk and sky-plane coordinates is given by: 

cos - sin \ / 1 \ / cos Oa -smOa \ / X 

sin 9i cos 9i cos^j — sin^j sin^a cos^a y 
1 / \ sinOi cos9i / \ 1 / \ z 

(4) 

Here, 9i is the angle in the sky plane from the X-axis to the ascending node of the disk on 
the sky (i.e. the point where a star orbiting in the disk plane has Z = and Z > 0). The 
inclination angle 9i is the angle between the normals to the sky plane and the disk plane 
(cos6'i = z ■ Z). The azimuthal angle 9a is measured in the disk plane from the ascending 
node of the disk on the sky to the symmetry axis of the disk (the positive x-axis) . All angles 
are measured in the usual sense, counterclockwise as seen from the positive Z, z, or z' axis. 

If the symmetry plane of the disk coincides with the symmetry plane of the large-scale 
M31 disk, then the parameters 9i and 9i are determined by the orientation of M31. The 
apparent major axis of the M31 disk has position angle PAgai = 37?7 (de Vaucouleurs 1958), 
and 9i — PAgai ± ^tt. The SE side of M31 is approaching, and the near side of M31 is to the 
NW (Hodge 1992), so the correct choice is 9i — PAgai — — — 52?3. The inchnation angle 
9i = 77?5 (Hodge 1992). 




2.2. Orbits 

In the models described here, we shall neglect the gravitational influence of the bulge 
and the nuclear disk on the disk-star orbits; that is, we assume that the disk stars travel 
on Keplerian orbits under the influence of a central BH of mass M,. This approximation is 
discussed more fully in §5. 
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The orbital period of a star traveling in a Keplerian orbit around the BH is P = 27r/n, 
where n — (///a^)^/^, fj, = GM, and a is the semimajor axis. Its coordinates {x',y',z' — 0) 
are given parametrically by the relations 

x' — a{cos E — e) , y' — aVl — sin£', M — E — esinE, (5) 

where e is the eccentricity, E is the eccentric anomaly, and M is the mean anomaly, defined 
as 2Trt/P where t is the time since periapsis passage. 

The velocity components of the star are given by 

. sinE . Vl - e2 cos E 

x' — —na — , y' — na — — , z — i). (bj 

1 — e cos E 1 — e cos E 

The corresponding linc-of-sight velocity V\os is then obtained by differentiating equations (3) 
and (4) with respect to time. 

We shall also use the canonical Delaunay variables (e.g. Dermott & Murray 1999), 

91 = M, Jl = y/JId, (7) 

92 ^ CO, J2 = y/iia{l - e2), (8) 
63 = Q, J3 = V/ia(l -e2)cos7. (9) 

The longitude of periapsis is zu = fl+uj. It is also convenient to introduce the eccentricity 
vector e, with components in the {x,y,z) coordinate system (ecosro, esinro, 0). 



2.3. The distribution function 

The nuclear disk is described by its distribution function /(x, v) , defined so that /(x, v)dx dv 
is the mass or light contained in the phase-space volume element dx dv. The corresponding 
volume element in orbital elements is dOdJ = ^fi^^'^a^^'^esmldadedl. According to Jeans's 
theorem, the distribution function (hereafter DF) can only depend on the integrals of motion, 
which in a Kepler potential can be taken to be a. e, /, Q. Since the DF is a function of five 
variables (only two of the components of e arc independent, since e lies in the z = plane), 
and we can only measure three (the sky-plane coordinates X and Y, and the line-of-sight 
velocity Vios), we cannot, even in principle, completely determine the DF directly from the 
observations. Instead, we choose a plausible parametric form for the DF, whose parameters 
are then fit to the observations. An alternative non-parametric approach would be to find 
the maximum-entropy DF that is consistent with the observations. 
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Thus we assume that 



f{a,e,I) = c/(a)exp 



em(a)] 



P 1 



(10) 



exp 



2(7/(a)2 



where em(a) = em(a)x is the mean eccentricity vector (note that e„j can have either sign). 
The magnitude of the mean eccentricity vector depends on semimajor axis but its direction 
is assumed to be fixed; this is consistent with the observation that gas-free stellar disks 
generally do not exhibit spiral structure, and with the anti-spiral theorem (Lynden-Bell & 
Ostriker 1967). Note that is the dispersion in one component of the eccentricity vector, 
not the dispersion in eccentricity; thiis, for = 0, the rms eccentricity is V2ae- Similarly, 
for small inclinations the rms inclination is \/2ai. For small eccentricity and inclination, the 
rms thickness at a given radius is simply (^^)^/^ = aai{a). 

For simplicity, the width of the eccentricity distribution, a^, is assumed to be indepen- 
dent of semimajor axis. However, if the disk is aligned with the plane of M31 and the width 
of the inclination distribution, aj{a), is independent of semimajor axis (so that the disk 
thickness is proportional to radius), the models are too flattened near P2. Thus we choose 
the form 



This variation of ai{a) with semimajor axis a is much less important for models that are 
not aligned with the plane of M31. In such models we have found that the best-fit value of 
aj tends to be large enough that the dependence of (J/(a) on semimajor axis is unimportant 
within the nucleus. 

Our numerical experiments show that the dispersion in inclination cTj has a profound 
effect on the photometric appearance of the disk at radii < 1". First consider aligned models. 
If the disk is too thin, the isophotes are too fiattened. If the disk is too thick, the surface 
brightness at PI is too low (i.e. smaller than at P2), and the "valley" between the peaks in 
the surface brightness distribution gets filled up. In the non-aligned model, if the disk is too 
thin, the edges of PI appear too well-defined, and the PI contours are crescent-shaped, in 
contrast to the more diffuse blob-like appearance of PI in the HST photometry; also, there 
is essentially no P2 peak, again in contradiction to the data. If the disk is too thick, the 
effect on the photometry of a non-aligned disk is the same as for an aligned disk. The effect 
of (jj on the kinematics is less dramatic, but still important: in the aligned model, for larger 
values of a^, the peaks of the rotation curve are further apart, the maximum rotation speed 
is higher, the rotation curve is more symmetric, and the displacement of the dispersion peak 
from P2B is smaller. In the non-aligned model, if the disk is too thin, the velocity extrema 
are too large compared to the data, and the dispersion peak is too narrow; if it is too thick, 
the converse happens. 



<7j{a) = CTj exp (—a/af) . 



(11) 



- 12 - 



We experimented with many forms for the mean eccentricity distribution 6^(0), includ- 
ing Gaussian, hnear, and exponential forms. If the eccentricity peaks at the origin, we found 
that the valley in the surface brightness profile between PI and P2 is not deep enough, so 
that PI and P2 do not appear as distinct as they do in the data. Thus, we concluded that the 
mean eccentricity should peak away from the origin. As suggested in the discussion following 
equation (2), we also found that satisfactory models required that the eccentricity decline 
rapidly to near zero in the semimajor axis interval from about 0'/5 to I'/O. We eventually 
chose the functional form 

[a-agf 



^mio) — Oi{ae — a) exp 



2^2 



(12) 



The functional form (10) that we have chosen for the DF allows eccentricities that exceed 
unity, which are unphysical; thus, in choosing orbital elements in our Monte Carlo simulation 
we discard points with e > 1. Figure 1 shows for the best-fit aligned and non-aligned 
models, along with the l-cx deviations represented by a^- A similar figure is shown by BOl 
(Fig. 22); they find a similar peak eccentricity ~ 0.7, although at a slightly larger radius 
(~ 0''5 instead of 0''3), and a similar sharp decline in eccentricity outside the peak. 

To choose the function g{a) that specifies the semimajor axis distribution, we first 
convert the DF (10) into an effective surface density distribution; here the effective surface 
density E(a) is defined so that 27raS(a)(ia is the mass with semimajor axes in the range 
[a,a-\-da\. We have 

E(ao) = / dJdeS{a-ao)f{J,e), (13) 



27ran 



where the action-angle variables (J, 6) are defined in equation (9). For the DF (10) we have 



S(ao) 



3/2 



1/2 



g{ao) / dl sin /exp 



2ai(aoy 



2n 



dw I dee exp 
'0 



[e - em{ao)Y 

2(7? 



(14) 

The innermost two integrals are easily evaluated by writing e de dw = de and using Cartesian 
coordinates for e: 



S(ao) 



1/2 



9(ao) [ 
Jo 



dl sin I exp 



2a,(ao)2 



(15) 



Since (7/(ao) is generally small, to a good approximation we have 



S(ao) = 271"^ ii^/'^ala]{ao)aQ^^'^ g{ao 



(16) 



Initially, we tried a simple exponential for the effective surface-density distribution, but 
we found that this form results in too much light at the center, and too much contamination 
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Fig. 1. — The mean eccentricity e^(a) as a function of semimajor axis, for the ahgned (left) 

and non-aUgncd (right) models. The gray bands represent 1-a deviations (parameter cTe). 
Note that Monte-Carlo points with |e| > 1 are discarded. The dashed line is the locus of 
orbits with apoapsis equal to O'.'S, the separation of PI and P2. 
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of the bulge by the disk at large radii. After some experimentation, we concluded that to 
match the photometry we require an effective surface-density distribution that has a central 
minimum, an outer cut-off, and a maximum roughly at the PI position. The steepness of 
the inner gradient of the distribution is important. For example, we experimented with a 
Gaussian distribution, S(a) = Sq exp[— A;(a — Oq)^], and an exponential model with a central 
hole, S(a) = So[exp(— a/ao) — exp(— a/ai)]. We found that for the former distribution the 
best-fit model (viewed from above the disk) has a wider central hole with a more well- 
defined boundary, a more well-defined PI peak which extends for a greater distance around 
the circumference of the disk, and less surface brightness coming from the outer extremities 
of the disk. When projected onto the sky plane, the Gaussian model is less successful in 
reproducing the dip in the surface brightness profile between the PI and P2 peaks. The 
fact that the double-exponential model is asymmetric about its maximum also seems to 
help in reproducing the photometric data. In the case of the double-exponential model, we 
consistently found best-fit models where ai ~ aoi leading to an inner gradient rising as a. 
We also experimented with models with S(a) = Soa" exp(— a/ao), and settled on n = 2 as 
the most promising value of this parameter. Therefore, we choose the following distribution, 
with a hole in the middle and an outer Fermi-type cutoff, 

1 exp [ci(a - C2)J 

the fit is quite insensitive to the value of the cut-off parameter ci. We then define g{a) using 
equation (16). 

In total, our model for the DF depends on 10 parameters: two for the semimajor 
axis distribution (ao, ci, C2); two for the inclination distribution ((7°, a/), and five for the 
eccentricity distribution (a, a^, Og, CTg). 



2.4. The bulge 

The M31 bulge is assumed to be spherical, isotropic and non-rotating. This relatively 
crude model is adequate because the bulge contributes only a modest background correction 

in the region R < 2" that is dominated by the nucleus. Wc use the photometric parameters 
of KB99 for our bulge model. At radii < 300", they found that the surface brightness of 
the bulge is well-described by a Sersic (1968) law, I{R) = /q cxp[— (i^/i^n)^/"], with Iq = 
15.40 mag arcsec"^ in V, Rn = 14 '.'0, and n = 2.19. With this bulge-nucleus decomposition, 
the central surface brightness of the bulge is 16% of the peak surface brightness of the 
nucleus. 
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Much of our analysis will be based on the bulge-subtracted kinematic data of KB99, in 
which case the bulge kinematics are irrelevant for our models. However, for some purposes 
we compare the model to kinematic data that includes the bulge light, and in this case we 
must model the bulge kinematics. To do so, we assume that the bulge has constant mass- 
to-light ratio T, compute the gravitational potential of the spherical Sersic profile described 
above, and add a central point mass corresponding to the (unknown) total mass of the BH 
plus nuclear disk. We then solve for the line-of-sight velocity distribution (LOSVD) at each 
radius, using formulae from Binney & Tremaine (1987) and Simonneau & Prada (1999). 

Wc assume that the mass-to-light ratios T of the bulge and nuclear disk are the same; 
this assumption is unlikely to be completely accurate since the colors of the nucleus and 
bulge are different, but the model predictions are only weakly dependent on the mass-to- 
light ratio of the bulge, and the mass-to-hght ratio of the disk only affects the model through 
the displacement between the BH and the bulge center, which is at the center of mass of the 
BH and the stellar disk. 



Numerical methods 



3.1. Constructing the nuclear disk 



We simulate the nuclear disk using a Monte-Carlo method, typically with = 10^ 
particles (there is little point in larger simulations, since is already comparable to the 
number of stars in the actual disk, and much larger than the number of giant stars, which 
dominate the light). First we assign semimajor axes a to the particles so that the distribution 
is uniform in T,(a')a' da' . For a given semimajor axis, the DF (10) implies that the mass 
in a small interval of the other orbital elements is 



dm oc exp 



[e 



(«)] 
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exp 



2ai{a) 



desinl dl dM dO.. 



(18) 



Thus, for each particle we choose Q and M from uniform probability distributions over the 
range [0, 27r); Cx and Cy from Gaussian distributions with means em(a) and respectively, and 
standard deviation (7^; and / from a Gaussian distribution weighted by sin/, with standard 
deviation ai{a). We obtain the eccentric anomaly E by numerical solution of Kepler's 
equation [the third of equations (5)]. We then use equations (5) and (6) to determine the 
position {x', y', z') and velocity {x', y', z'), and equations (3) and (4) and their time derivatives 
to determine the positions and velocities in the sky-plane coordinate system. 
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3.2. Photometry 

We fit the model to the deconvolved HST V-band (F555W) image presented in L98. 
The intensity scale is calibrated using the zero-point given in L98. We assume that the pixel 
containing the maximum luminosity in the blue cluster in the P2 peak is the position of 
the BH, and this point is also the origin of our coordinate systems. The excess brightness 
contributed by the UV cluster at P2 in the data is clipped out prior to fitting, since we may 
assume that these stars have a small mass-to-light ratio. We compare the data to the model 
on a sky-plane grid with cell size equal to the pixel size, 0'/0228 = 0.085 pc, and a circular 
boundary of radius Rg — 5 pc — 1'.'34. 

We add the surface-brightness contribution from the bulge by assuming that the center 
of the bulge coincides with the center of mass of the disk-BH system; thus the center of the 
bulge is not precisely at the origin. 



3.3. Kinematics 

We fit the kinematics to the bulge-subtracted spectroscopic data from KB99, which 
combine high S/N with reasonably good spatial resolution. These data have a scale of 
0'/0864/pixel, a PSF with FWHM of 0'/64, and a slit width 2s = 0'/35. Two spectra were 
taken at position angle PA = 50° and two more at PA = 55°. KB99 then shifted these four 
spectra to a common center and co-added them for further analysis. 

We reproduce this procedure in the model. We compute the LOSVD from the nucleus 
on a sky-plane grid with cell size equal to the pixel size, and then smooth the LOSVDs by 
convolving with the PSF given in equation (3) of KB99 and a top-hat slit of width 2s. We 
then add the LOSVDs at PA = 50° and PA = 55°. Following KB99, we fit the LOSVDs 
with a fourth-order Gauss- Hermite expansion (van der Marel & Franx 1993; Gerhard 1993): 



(19) 



The coefficients and parametrize the lowest order odd and even deviations from Gaus- 
sian line profiles. The velocity bins are 5 km s^^ wide and are weighted equally in the fit. 
The kinematic parameters extracted from the fit (V, a, /13, /14) can then be compared with 
the observations from KB99.^ 



•^Note that the kinematic data in KB99 (their Table 2) have as their origin the position where the rotation 
speed V = relative to the systemic velocity of M31, which is shifted from our origin at the center of P2B 
toward PI by 0'.'051 ± 0.014 according to KB99 or O'.'OSl according to BOl. 
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3.4. Pcirameter fitting 

We fit the sky-plane image to the model using the statistic 

Xi = -^^ [datagrid(i, j) - A modelgrid(i, j)f , (20) 

where A specifies the overall normahzation of the model, w{i,j) is the weight assigned to 
pixel (i, j), and Rg is the radius of the region we are fitting. 

One potential drawback of this approach is that most of the weight comes from the outer 
parts of the image, while the most critical photometric test of the model is its ability to fit 
the observations near the center, in particular along the P1-P2 axis. We have addressed 
this problem in two ways: (i) we have used a weight function w{i,j) — •\/datagrid(i, j), 
which gives more weight to the high-surface brightness regions in PI and P2 (the form of 
this weighting factor was selected after experimenting with a variety of possibilities), (ii) 
We have computed a second statistic, xl; which is similar to but restricted to pixels in 
a slit (of width one pixel on the model photometric grid) that extends through PI and P2 
to a distance ±-Rg from the origin. This statistic focuses on the difficult task of fitting the 
double structure of the nucleus. 

We specify the fit to the kinematic data by a statistic Xg, which measures the mean- 
square differences of V and a between the model and the KB99 bulge-subtracted data at 
46 points along the PA=50°-55°axis (Table 3 of KB99). Each data point is weighted by the 
observational uncertainty given by KB99. Note that KB99 present velocity and dispersion 
measurements extracted with two techniques: Fourier Quotient (FQ, their Tables 2 and 
3) and Fourier Correlation Quotient (FCQ, their Tables 4 and 5). The Gauss-Hermite 
fitting procedure outhned above mimics the FCQ algorithm more closely; however, the FQ 
data are given at significantly finer spatial resolution. The velocity differences between the 
two algorithms arc < 2 km s~^; the close agreement between the two sets of velocity and 
dispersion data is illustrated in Figures 5 and 6. Thus we choose to fit to the KB99 FQ data. 

We determine the best-fit model by minimizing the statistic -|- A1X2 + ^^xi^ where 
the relative weights Ai and A2 are chosen by hand such that for the initial state model for 
the full photometry and kinematics fit, Xi ~ -^1X2 ^"^^ (Xi + ^1X2) ~ ^sxi- The procedure 
used to determine the initial state model is described later in this section. 

We do not fit to KB99's measurements of the Gauss-Hermite parameters and h^, 
or to spectroscopic data from other investigators. However, we shall compare the best-fit 
models we obtain in this way to STIS observations of kinematic parameters including and 

/l4. 
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To reduce the dimensionality of our large parameter space, we choose two parameters 
by hand: Ci and T = M/Ly. The Fermi cut-off parameter ci (eq. 17) does not significantly 
affect the fit within the inner 2", and the photometry is quite insensitive to its exact value, 
so we simply choose Ci = 4 pc~^. The mass-to-light ratio T does not enter the kinematic fit 
since we use bulge-subtracted kinematic data, and only affects the photometric fit through 
the displacement between the BH and the bulge center, which is the barycenter of the BH 
and nuclear disk. To find the position of the barycenter, we simply assume the same mass- 
to-hght ratio for the nuclear stars as the KB99 value for the bulge stars, Ty = 5.7. We 
then use the resulting disk mass and the model's BH mass to compute the position of the 
barycenter. 

If the nuclear disk is assumed to be ahgned with the M31 disk ( "aligned models" ) we 
now have 11 free parameters: one orientation parameter {Oa, the angle in the disk plane 
between the sky plane and the symmetry axis of the nuclear disk); one mass parameter 
(the BH mass M,), and 9 parameters for the disk DF (oq, C2, ctj, a/, a, Og, Oe, w, (Je). If the 
nuclear disk orientation is fitted to the data ("non-aligned models"), we must fit two more 
orientation parameters, Oi and Oi. 

The fitting is done with a downhill simplex algorithm (Press et al. 1992). 

To obtain initial conditions for the fitting program, we made use of the thin disk model. 
It is relatively easy to find parameters for equations (12) and (17) that give a thin disk 
X-axis surface-brightness profile (1) similar to the observed profile along the P1-P2 axis. 
Hence, reasoning that thickening the disk would lower the surface brightness at PI, we 
chose as a starting point an "extreme" thin disk model which maximized the ratio of surface 
brightnesses P1/P2. 

This extreme thin disk model produced optimum values for the parameters [oq, Oi, a, a^,, Og, w] . 
These were then used as a first guess to fit to the photometry of the three-dimensional disk 
model (using only the statistic xi)-, with plausible initial guesses for aj, 9a, cr'j and o"e- When 
the fitting algorithm converged for this limited parameter set, they were then input into the 
full kinematic and photometric fit as an initial guess, adding M, to the fit. 

It should be emphasized that the thin disk model was only used to obtain the initial 
conditions; the rest of the fitting process uses the thick disk model of §2. 



X (arcsec) X (arcsec) 

Fig. 2. — Surface-brightness distribution in the nucleus of M31. The left panel shows the 
data (V-band; from L98) and the right panel shows the best aligned model. Contours are at 
0.25 mag intervals. North is 55.7° to the left of the top of the plot. 




Fig. 3. — Surface-brightness distribution in the nucleus of M31. The left panel shows the 
data (V-band; from L98) and the right panel shows the best non-aligned model. Contours 
are at 0.25 mag intervals. North is 55.7° to the left of the top of the plot. 
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Fig. 4. — Surface brightness along a slit of width 1 pixel (= 0'/0228) that extends through 
PI and P2. The solid line represents the data (from L98) and the dashed lines are the best 
aligned model (left panel) and non-aligned model (right panel). 




r (arcsec) r (arcsec) 

Fig. 5. — Rotation speed V along the axis at PA=50°-55°, as determined by fitting the 
LOSVD to the Gauss-Hermite expansion (19). The data from KB99 after bulge subtraction 
are shown as black dots (FCQ) and crosses (FQ), and the solid lines are the best aligned 
model (left panel) and non-aligned model (right panel). 
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Fig. 6. — The dispersion a along the axis at PA=50°-55°, as determined by fitting the 
LOSVD to the Gauss-Hermite expansion (19). The data from KB99 after bulge subtraction 
are shown as black dots (FCQ) and crosses (FQ), and the solid lines are the best aligned 
model (left panel) and non-aligned model (right panel). 




Fig. 7. — The Gauss-Hcrmitc parameter h-^ along the axis at PA=50°-55°. The data from 
KB99 after bulge subtraction are shown as black dots (FCQ), and the solid lines are the best 
aligned model (left panel) and non-aligned model (right panel). 



-22- 



Table 1. Model parameters 



parameter 


aligned model 


non-aligned model 


a 


0.288 


0.197 


tte (pc) 


3.97 


4.45 


ag (pc) 


1.51 


1.71 


w (pc) 


1.53 


1.52 


M.(Mo) 


9.63 X 10^ 


1.02 X 10^ 


ao (pc) 


4.61 


1.37 


ci (pc-^) 


4a 


4a 


C2 (pc) 


3.79 


4.24 


O-e 


0.344 


0.307 


Oa (deg) 


-11.0 


-34.5 


a/ (pc) 


13.7 


31.5 


a° (deg) 


36.2 


24.6 


(deg) 


77.5^ 


54.1 


^« (deg) 


-52.3^ 


-42.8 


T(Mo/Lo) 


5.7^ 


5.7^ 



*This parameter was not fitted. 
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Fig. 8. — The Gauss-Hermite parameter /t4 along the axis at PA=50°-55°. The data from 
KB99 after bulge subtraction are shown as black dots (FCQ), and the sohd hues are the best 
aligned model (left panel) and non-aligned model (right panel). 



Table 2. Model photometry 



parameter 


A 


N 


data 


Separation of PI and P2 (arcsec) 


0.6 


0.5 


0.5 (L98) 


P1-P2 position angle (deg) 


39° 


42° 


43° 


Separation of P2 from bulge center (arcsec) 


0.04 


0.04 


0.07 (KB99) 


Disk luminosity < 1" projected {Lq) 


2.9 X 10^ 


2.9 X 10^ 


2.9 X 10^ 


Disk luminosity inside 4" x 4" box {Lq) 


3.6 X 10^ 


3.6 X 10^ 


5.4 X 10^ 



Note. — A=aligned model; N=non-aligned model. Disk luminosity is computed 
using bulge model from KB99. The projected model disk luminosity is normalized 
to the projected bulge-subtracted data luminosity < 1". 
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4. Results 

We shall describe two models: the best-fit aligned model, in which the nuclear disk is 
assumed to lie in the symmetry plane of the M31 disk, and the best-fit non-aligned model, 
in which the orientation of the nuclear disk is chosen to optimize the fit. 

The aligned model we describe here was obtained by minimizing the statistic Xi/(6 x 
10^) + x\/{^ ^ 10^) + xi) while the non-aligned model was based on the statistic x^/(1.2 x 
10'') -|- x|/(1.2 X 10^) + x\- The parameters of the models are given in Table 1. 

Figures 2 and 3 show the surface-brightness distributions in the data and the two mod- 
els. Both models correctly reproduce the apparent double structure of the nucleus, and the 
approximate sizes and shapes of PI and P2. However, in the aligned model the maximum 
surface brightness of PI is lower than the maximum brightness in P2, while in the observa- 
tions PI is brighter. The non-aligned model does much better in reproducing the relative 
surface brightnesses of PI and P2. 

The position angle of the line joining PI and P2 is 39° in the aligned model and 42° in 
the non-aligned model, close to the observed value of 43°. The position angle of the major 
axis of PI in the data (as determined by the isophotes with semimajor axis less than 0'/2) 
is 63°. This 20° isophote twist is not reproduced in the models. However, the observed 
twist is strongly affected by the limited number of giant stars that contribute to the surface 
brightness in the middle of PI. We have explored this effect by reducing the number of stars 
in our Monte-Carlo realization by a factor ~ 100 to approximate the actual number of giants 
in the nucleus (of order 10^). We find that in this case the position angle of the PI major 
axis varies significantly in different realizations; therefore, the difference in position angle 
between the major axis of PI and the P1-P2 axis does not provide much useful information. 
However, the range of position angles seen in realizations of the non-aligned model is larger 
than in the aligned model, so that the non-aligned model is more easily able to accommodate 
the 20° isophote twist seen in the data. 

Figure 4 shows the surface brightness along the P1-P2 axis in both models. In the data, 
the surface-brightness peak at PI is approximately 0.4 mag brighter than the peak at P2. 
In the aligned model, PI is fainter than P2 by 0.1 mag, while in the non-ahgned model PI is 
brighter by 0.2 mag; thus the fit of the non-ahgned model to the relative brightnesses of PI 
and P2, though not perfect, is much better than the aligned model. The brightness of PI 
can be increased in the aligned model by reducing the disk thickness (smaller ctj) but then 
PI becomes too elongated. 

The non-ahgned model correctly reproduces the dip in surface brightness between PI 
and P2, but the minimum is displaced 0'/2 further toward P2 than it is in the data. A feature 
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of both models is that the P2 peak is displaced to the anti-Pl side of the our origin at the 
BH (i.e. the UV peak P2B). This is a consequence of the luminosity depression imposed at 
the origin (S(a) oc a?, eq. 17). While in principle this displacement could be an important 
test of the model, with implications for the kinematics (specifically, the displacement of 
the dispersion peak to the anti-Pl side), we have found that this displacement cannot be 
determined reliably if the number of stars in the simulation is reduced to approximate the 
number of giant stars in the actual system: in this case Monte-Carlo simulations can be 
created which have no significant displacement of the P2 peak from P2B. 

Note that both models fail to reproduce the surface brightness outside ~ 1'.'3, because 
the surface brightness in the model falls off more sharply than in the data. The obvious 
fix is to increase the cutoff parameter C2 in equation (17). However, this change affects the 
fit to the kinematics, decreasing the amplitude of the rotation minimum and lowering the 
dispersion peak. The discrepancy in the surface brightness at large radii could arise either 
because the real bulge is cuspier than our model, which we adopted from KB99, or because 
the nucleus contains extended wings that are not captured by our parametrization. 

Figure 5 shows the rotation speed V — more precisely, the parameter V in the Gauss- 
Hermite expansion (19) — along the average of the two shts at PA=50° and 55° measured in 
KB99. The aligned model, shown in the left panel, reproduces the approximate shape of the 
rotation curve, in particular the strong asymmetry (the peak of the observed rotation curve 
is at 179 km s^^ on the PI side and 236km s^^ on the anti-Pl side; see Table 3 for the peak 
velocities in the models). The zero-point of the model rotation curve is displaced from P2B 
toward PI by 0'.'05, consistent with the displacement of O'.'05-O'.'IO estimated by KB99. The 
most significant discrepancy is that the maximum rotation speed on the anti-Pl side, near 
0'/6, is too small by about 24 km s"^ or 10%; the rotation speed near the peak on the PI 
side is also too low, but by a smaller amount. 

The right panel shows the rotation curve for the non-aligned model, which fits the 
rotation curve substantially better than the best aligned model. 

Figure 6 shows the dispersion profile along the same slits. The dispersion profile in the 
aligned model, shown on the left, fits the data fairly well, but (i) the maximum dispersion is 
too low by 14% (248km s""*^ compared to 287km s"^); (n) the model dispersion is too large 
by about 10km s~^ on the PI side between radii of 0'/4 and I'.'l. The dispersion is also too 
high near 2" on the anti-Pl side, but this discrepancy probably is due to the shortcomings 
of the bulge model discussed above. 

The model correctly reproduces the displacement of the dispersion peak from P2B by 
about 0'/13 in the anti-Pl direction — the displacement in the model is about 0'/17. We could 
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increase the peak dispersion in the model by increasing the BH mass, but then the excess 
model dispersion on the PI side would become even worse. 

We have attempted to increase the maximum dispersion by adding a compact spherical 
stellar system at P2 (in the form of a Plummer sphere) ; the hope was that the high- velocity 
stars in a small, low-luminosity system of this kind would increase the dispersion without 
degrading the fit to the photometry. However, we found that increasing the dispersion profile 
in this manner made the rotation curve too symmetrical, so the overall fit was not improved. 

The dispersion profile and rotation curve of the aligned model could be made to fit the 
KB99 data if we assume that the width of the PSF is 10-15% smaller than the value used 
by KB99. However, it is unlikely that there is an error of this magnitude in the estimated 
PSF. 

The right panel of Figure 6 shows the dispersion profile in the non-aligned model. The 
model profile fits the data substantially better than the best aligned model, though the 
maximum dispersion is still too low by 5%. The dispersion peak is displaced by O'.'OQ from 
P2B in the anti-Pl direction, in good agreement with the observed displacement of 0'.'13. 

Figures 7 and 8 show the parameters hs and /i4 in the Gauss-Hermite expansion. Both 
the aligned and non-aligned models are moderately successful at reproducing the data within 
the large uncertainties; thus it appears that a model that is fit to the overall scale and center 
of the velocity distribution (measured by a and V) also reproduces the gross features of 
its shape. Note that we do not use data on cither h^, or /?.4 in the fitting process, so the 
comparisons in these figures are a measure of the predictive power of the eccentric-disk 
model rather than the quality of the fit. 

BOl have used the OASIS integral-field spectrograph, assisted by adaptive optics, to 
measure the two-dimensional mean- velocity and velocity-dispersion fields in the central I'.'S 
of M31, with a PSF having FWHM of ~0'.'4-0:'5. They find that the kinematic axis (the 
axis of symmetry of the mean- velocity field) has a position angle of 56°, which is rotated 
from the P1-P2 axis at position angle 43°. Figures 9 and 10 show the mean- velocity and 
velocity-dispersion fields for the ahgned and non-aligned models, as viewed with a Gaussian 
PSF having FWHM of 0'.'13 to smooth out the shot noise in the contours. These figures 
can be directly compared to the model photometry in Figures 2 and 3, which has the same 
resolution and orientation. Once again, the non-aligned model provides a better fit to the 
data than the aligned model. In particular, the kinematic major axis (defined by the line 
joining the velocity extrema) in the non-aligned model has PA=48-50°, much closer to the 
observed orientation of 56° than the kinematic major axis in the aligned model, PA=40°. 
We suspect that the orientation of the kinematic axis is resolution-dependent, which may 
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Table 3. Model kinematics 



parameter 


A,S 


N,S 


STIS (FCQ) 


A,K 


N,K 


KB99 (FQ) 


PI rotation maximum (km s""*^) 


193 


198 


203 


170 


178 


179 


Position of rotation maximum (arcsec) 


-0.76 


-0.61 


-0.59 


-0.95 


-0.86 


-0.99 


Anti-Pl rotation minimum (km s~^) 


-277 


-369 


-319 


-212 


-228 


-236 


Position of rotation minimum (arcsec) 


0.41 


0.20 


0.25 


0.78 


0.60 


0.65 


V — position (arcsec) 


-0.17 


-0.12 


-0.12 


-0.05 


-0.05 


-0.05 


Dispersion peak (km s~^) 


252 


328 


302 


248 


274 


287 


Position of dispersion peak (arcsec) 


0.15 


0.10 


0.08 


0.17 


0.09 


0.13 



Note. — A=aligned model; N=non-aligned model; K=KB99 resolution; S=STIS resolution 




X (orcsec) X (orcsec) 

Fig. 9. — Two-dimensional kinematics of the best aligned model. The left panel shows the 
mean-velocity contours and the right panel shows the velocity-dispersion contours; both pan- 
els are superimposed on a gray-scale image of the nucleus. The contour interval is 25 km s~^. 
The thick hue is the zero (systemic) velocity contour (left panel) and the 200 km disper- 
sion contour (right panel). The locations of the center of PI (filled circle) and P2B (filled 
diamond at the origin) are shown. North is 55.7° to the left of the top of the plot. 
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account for the difference of 7° between the non-ahgned model and the BOl data. In Figure 
7 of BOl, the isovelocity contours between PI and P2 are nearly perpendicular to the line 
joining the velocity extrema, while this angle is oblique in our non-aligned model. We find 
that lower-resolution maps of the velocity field of the non-aligned model yield isovelocity 
contours that are more nearly perpendicular to the kinematic axis, consistent with the BOl 
data. 



4.1. Comparison with STIS data 

We now compare the predictions of the models described above to new high-resolution 
STIS observations (Bender et al. 2003). Recall that our models were fitted to the ground- 
based spectroscopic data from KB99; here, the kinematics have simply been recomputed 
for STIS resolution (slit width O'.'l) without modifying the fit. The effects of the instru- 
mental broadening on the LOSVD have been accounted for using a template star observed 
through the same slit (R. Bender, private communication). We use the PSF given by Bower 
et al. (2001), which can be modeled as the sum of two Gaussians (K. Gebhardt, private 
communication) : 



G{r) = 1.01 exp 



r.2 



„ (r- 0.105)^ 



2(0.0329)2 



+ 0.12 exp 



2(0.0384)2 



(21) 



where r is the distance in arcseconds. Note that a Wiener filter has been applied to the STIS 
data (R. Bender, private communication) to remove high frequencies, a step that we have 
not taken in analyzing the model kinematics. 

The spatial resolution of the STIS observations is far better than that of the KB99 
observations (the slit width is a factor of 3.5 smaller, and the FWHM of the PSF is smaller 
by a factor of 5.4). Thus the comparison of the STIS kinematics with our models is a 
stringent test of their predictive power. 

The kinematics of the best-fit aligned and non-aligned models that would be observed 
with STIS resolution are shown in Figures 11-14. The solid curves are bulge-subtracted and 
the dashed curves include the bulge. The data, shown as dots (FCQ) and crosses (FQ), are 
obtained at PA=39° and are bulge-subtracted. 

In contrast to the KB99 data, the FQ and FCQ results are substantially different for 
the STIS data. The difference arises because FQ assumes Gaussian LOSVDs while FCQ 
simultaneously derives velocities, dispersions, hs, and by fitting the LOSVD with a Gauss- 
Hermite series. Thus, the results from the two techniques do not agree where LOSVDs have 
strong wings and the higher order Gauss-Hermite moments are large. Since we mimic the 
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FCQ procedure in our calculations, one should compare the model curve with the FCQ 
results (the dots). 

The non-aligned model is remarkably successful at reproducing the rotation curve on 
the PI side, as well as the rotation gradient through the center (right panel of Fig. 11). On 
the anti-Pl side the peak rotation speed in the model is too high by 16% (see Table 3) and 
the model rotation curve outside the peak is also too high. 

The displacement of the zero of the rotation curve is 0'.'12 toward PI in both the data 
and the non-aligned model. 

The dispersion peak in the non-aligned model is 9% too high (Fig. 12). The dispersion 
peak position is reproduced correctly — in the data the peak is displaced in the anti-Pl 
direction by O'.'OS, compared to O'/IO in the model. 

At this resolution, both models show a secondary peak in the dispersion curve (Fig. 
12). This peak is a consequence of the hole in the middle of the model surface density 
distribution. The data do not conclusively favor or disfavor a feature at this location. One can 
possibly make this feature less prominent by making the central surface-brightness gradient 
shallower. However, we have been unable to remove the central hole completely without 
drastically worsening the fit to the photometry. Thus, wc believe that a secondary peak in 
the dispersion profile should generally be present at some level in the type of model proposed 
here. 

The model dispersions are much lower than the dispersions in the data at radii around 
1". We believe that this discrepancy arises because the KB99 bulge model and our disk 
model do not have enough flexibility to match the photometry and kinematics in this region; 
however, this problem should not affect our models at smaller radii. 

Both the aligned and non-aligned models are quite successful in reproducing the overall 
shape of the profiles of the Gauss-Hermite parameters and /14 (Figs. 13 and 14). In 
particular, the models reproduce the shape, amplitude, and location of the sharp dip in 
/i3 and the peak in /i4. The agreement between data and models is particularly impressive 
since the model is derived by fitting only to the low-resolution mean velocity and dispersion 
profiles. 

Figure 15 shows the LOSVDs for each model at a few locations near P2. These LOSVDs 
are calculated at STIS resolution along PA=39°, and do not include the bulge. Near the 
BH, the LOSVDs are very asymmetric and have strong wings on the prograde side of the 
mean rotation velocity. These features are also seen in the STIS LOSVDs (plotted here in 
long dashes). However the model wings are more extensive than in the data. 
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Note that the model LOSVDs are constrained to be positive, while the STIS LOSVDs 
are not; thus the latter show (unphysical) low-amplitude oscillations that are not present in 
the models. 



4.2. Internal dynamics of the disk 

T95 suggested that two-body relaxation would lead to significant thickening of a thin 
disk if the disk age were comparable to the Hubble time. In this subsection we investigate 
this process quantitatively for our disk models. For simplicity, we shall neglect the mean 
eccentricity of the disk, that is, we approximate the disk as axisymmetric. 

Axisymmetric gravitational stability in a razor-thin disk requires that Toomre's Q pa- 
rameter exceeds unity, where (Toomre 1964; Binney & Tremaine 1987) 

Q = 0.30(7e-^; (22) 

here E(a) is the surface density at radius a and is the rms value of one component of 
the eccentricity vector or 2~^^'^ of the rms scalar eccentricity. Using the surface-density 
distributions in the aligned and non-aligned models, we find that the value of CTg required 
for stability never exceeds 0.1 (Fig. 16); once the non-zero disk thickness is included, the 
required value would be even lower. Thus the aligned and non-aligned disk models, which 
have (7e = 0.34 and 0.31, are safely stable. 

Let us assume that the disk formed to — 1 x 10^° yr ago, and was initially cold (i.e. (jg 
just large enough for gravitational stability). Two-body relaxation will gradually increase 
both the rms eccentricity and the rms inclination of the disk stars. This process of "disk 
heating" or "viscous stirring" has been studied thoroughly in protoplanetary disks. It is 
found that disk heating leads to a specific ratio between the rms inclination and eccentricity, 
ai ~ 0.5(Te, and that (Ohtsuki, Stewart, & Ida 2002) 

(O^(to) = 4.52Qto^^^lnA, A^O.G^JlY'; (23) 

where (P) — 2a] is the mean-square inclination, Jl(a) = {GM^/a^y^'^, and m* is the stellar 
mass. In deriving this formula we have neglected the time variation of (P) in the factor A, 
since it only appears in a logarithm, and we have assumed that the disk is in the dispersion- 
dominated regime, where {P)^/"^ ^ (m*/M,)^/^ ~ 0.002. 

The predicted eccentricity and inclination dispersions, ui — {PY^'^ / y/2 and Ue — 2a i, 
are shown in Figure 16 for the aligned and non-aligned models. We have shown both the 
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Fig. 10. — Two-dimensional kinematics of the best non-aligned model. Parameters as in 
Figure 9. 
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Fig. 11. — Rotation speed V along the axis at PA=39°, as observed with STIS resolution 
(see text). The rotation speed is determined by fitting the LOSVD to the Gauss- Hermite 
expansion (19). The best ahgned model is shown in the left panel and the best non-ahgned 
model is shown in the right panel. The solid curve is bulge-subtracted and the dashed curve 
includes the bulge. The bulge-subtracted STIS data from Bender et al. (2003) are shown as 
black dots (FCQ) and crosses (FQ). 
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Fig. 12. — Velocity dispersion a along the axis at PA=39°, as observed with STIS resolution 
(see text), determined by fitting the LOSVD to the Gauss- Hermite expansion (19). The best 
aligned model is shown in the left panel and the best non-aligned model is shown in the right 
panel. The sohd curve is bulge-subtracted and the dashed curve includes the bulge. The 
bulge-subtracted STIS data from Bender et al. (2003) are shown as black dots (FCQ) and 
crosses (FQ). 
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Fig. 13. — The Gauss-Hcrmitc parameter along the axis at PA=39°, as observed at STIS 
resolution. The best aligned model is shown in the left panel and the best non-aligned 
model is shown in the right panel. The solid curve is bulge-subtracted and the dashed curve 
includes the bulge. The bulge-subtracted STIS data from Bender et al. (2003) are shown as 
black dots (FCQ). 
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Fig. 14. — The Gauss-Hcrmitc parameter /14 along the axis at PA=39°, as observed at STIS 
resolution. The best aligned model is shown in the left panel and the best non-aligned 
model is shown in the right panel. The solid curve is bulge-subtracted and the dashed curve 
includes the bulge. The bulge-subtracted STIS data from Bender et al. (2003) are shown as 
black dots (FCQ). 
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Fig. 15. — The LOSVD distributions near P2 at STIS resolution for the ahgned (dotted) 
and non-ahgned (soUd) models. The LOSVDs from the STIS data are plotted in long dashes 
(in places the STIS LOSVD is negative and these sections are not shown). The curves are 
normalized to unit area. The velocity zero point is the systemic velocity of M31. Note that 
not all the LOSVD bins have the STIS LOSVDs overplotted, because they were provided 
only at specific radii, and we could only match some of these radii to a corresponding pixel 
on our LOSVD grids. 
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Fig. 16. — The inclination and eccentricity dispersions ai (solid lines) and Ue (dashed hnes) 
in the aligned model (left panel) and the non-aligned model (right panel). The quantities 
2a| and 2a1 are the mean-square inclination and eccentricity at a given semimajor axis in a 
disk with mean eccentricity = 0. Also shown are the dispersions expected from two-body 
relaxation if the disk age is 10^'^ yr, from equation (23) and the relation o"/ = O.Sce. The 
curve labeled (7e{Q = 1) is the minimum eccentricity dispersion required for axisymmetric 
gravitational stabihty (eq. 22). 
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theoretical predictions derived from equation (23) and the surface density and BH mass from 
Table 1, and our empirical fits to the eccentricity and inclination dispersions, also from Table 
1. Note that the theoretical predictions should be robust, since the predicted value of o"/ 
depends on only the fourth root of the age, surface density, stellar mass, etc. (eq. 23); the 
principal uncertainty is probably our assumption that the disk is axisymmetric. 

Two-body relaxation requires that the disk axis ratio in the radius range 0.5-1" is at 
least 0.13-0.17 in the aligned model, and 0.23-0.25 in the non-aligned model (these values 
are a/, which equals the ratio {z^Y^"^ /r in a circular disk). In the aligned model, the thickness 
required by relaxation is much smaller than the model thickness (o"/ = 0.58-0.55 in the same 
radius range), and thinner disks are strongly excluded by the photometry. Thus, if this model 
is correct, some process other than two-body relaxation must determine the thickness. The 
non-aligned model is thinner, ct/ = 0.38-0.41, and the predicted thickness due to relaxation 
is larger, so it is possible that in this case the thickness is largely determined by two-body 
relaxation. In aligned disks there is good agreement between the eccentricity dispersion (Xe 
due to relaxation and the best-fit value in the model, while in the non-aligned model the 
eccentricity dispersion due to relaxation actually exceed the best-fit value in the model. We 
suspect that these disagreements are consistent with the uncertainties in the model. 



5. Discussion 

The extremely short dynamical time (1.5 pc/200km s~^ ^ 7500 yr) of the M31 nucleus 
strongly suggests that it is in dynamical equilibrium, an assumption that we adopt through- 
out this paper. The M31 nucleus reveals by far the richest phenomenology of any equilibrium 
stellar system, and hence provides a unique challenge for stellar dynamics. 

We have fitted parametrized eccentric-disk models of coUisionless stellar systems to the 
photometry and kinematics of the nucleus. These models explain most of the features in 
the nucleus successfully. However, despite the large number (11-13) of free parameters and 
our efforts to use plausible fitting functions for the DF, there remain significant differences 
between our best-fit models and the observational data: for example, we are unable to fit 
precisely the separation and relative surface brightnesses of the PI and P2 components (Fig. 
4). The important question of whether these discrepancies arise from shortcomings of the 
eccentric-disk hypothesis, approximations in our model (e.g. the orbits are assumed to be 
Keplerian), the limited number of bright stars in the nucleus, or the limited flexibility of our 
parametrization remains to be determined. 

The best-fit black-hole mass is 1.0 x IO^Mq (Table 1); this estimate is hkely to be about 
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10%-20% too high because we have neglected the contribution of the nuclear disk stars to the 
total mass [the ratio of the disk mass inside 1" to the BH mass is ~ 0.15 (Table 2)]. Previous 
estimates of the BH mass include 3-7 x lO^M© (Dressier & Richstone 1988), lO^-^-lO^M© 
(Kormendy 1988), 4-5 x WMq (Richstone, Bower & Dressier 1990), 7 x WMq (Bacon et 
al. 1994), 7.5 x IO^Mq (T95), 7-10 x WMq (Emsellem & Combes 1997), (3 ±1.5) x lO^M© 
(KB99), and (7+3^) x lO^M© (BOl). The relatively smaU mass preferred by KB99 is found 
by a quite different method than the other estimates (they use the displacement between 
P2B and the center of the bulge, as opposed to the rotation curve and velocity dispersion 
profile), and may be subject to systematic errors despite the great care taken by KB99 to 
measure this small displacement accurately. For comparison, the correlation between BH 
mass and bulge dispersion observed in other nearby galaxies (Tremaine et al. 2002) predicts 
a mass of 6 x IO^Mq with an uncertainty of a factor of two. 

Although the models presented in this paper represent a substantial advance in accuracy 
and realism, they still have many shortcomings. 

• We assume that the stars in the nuclear disk orbit in the point-mass potential of the 
BH, that is, we neglect the gravitational forces from the disk and bulge. Neglecting 
the bulge is reasonably safe: within 1" the mass of the bulge is < IO^Mq, less than 
1% of the mass of the BH. Neglecting the disk is a more serious approximation: in our 
best- fit models, the disk mass within 1" is ~ 15% of the BH mass (Tables 1 and 2); 
thus the effects of the self-gravity of the disk are small but not neghgible. Including 
the gravitational potential of the disk will modify the kinematics of the disk orbits, 
and may permit new types of orbit family not present in the spherically symmetric 
BH potential (e.g. lens orbits and chaotic orbits; see Sridhar & Touma 1999 and 
Sambhus & Sridhar 2002). An additional consequence of including the disk potential 
is that the eccentric disk will precess, an effect not included in our model (Sambhus & 
Sridhar 2000, 2002). We believe that properly including these effects may significantly 
improve the kinematic fits. However, this improvement is unlikely to have much effect 
on the photometry. Thus, in particular, our conclusion that non-aligned models fit the 
photometry better than aligned models is unlikely to change. 

• The eccentric-disk model requires that the apsides of the disk stars precess uniformly so 
that the disk maintains its apsidal alignment; this requirement provides an additional 
strong constraint on the mass and eccentricity distribution in the model that we have 
not exploited. Statler (1999) has argued that uniform precession requires that the 
disk have a steeply declining eccentricity gradient and a change in direction of the 
eccentricity vector. The steep eccentricity gradient is present in our models, but is 
required by the high surface brightness at PI rather than for any dynamical reason 



-39- 



(eq. 2). The eccentricity does not reverse sign to any significant extent (see Fig. 1). 
Statler's argument for the eccentricity reversal is plausible but suspect, for at least two 
reasons: (i) The same argument, applied to a slightly eccentric test-particle orbit in 
an axisymmetric disk, would predict that the periapsis would precess in a prograde 
direction, but orbits in continuous disks precess in a retrograde direction {w = Q — k, 
where fl and k are the azimuthal and radial frequencies, and in general k > fl). (ii) 
Tremaine (2001) has shown that self-consistent m — 1 hnear normal modes in nearly 
Keplerian disks require non-zero velocity dispersion or thickness, effects that are not 
included in Statler's numerical calculations or qualitative discussion. 

• Our models contain only prograde orbits, and retrograde orbits may also contribute to 
the structure of the M31 nucleus; Sambhus & Sridhar (2002) find that retrograde orbits 
significantly improve their fits and this possibility deserves to be explored further. 

• We have used parametrized models for the DF; thus, we are unable to determine 
whether the differences between our models and the data are due to shortcomings 
of the model or the parametrization. A more fiexible, but computationally expensive, 
approach would be to compile an orbit library and use a penalized maximum-likelihood 
method such as maximum entropy to construct the best-fit smooth DF. 

• The bulge model of KB99 does not match on well to our disk models, and there are 
discrepancies in the photometry (Fig. 4) and the velocity dispersion (Figs. 6 and 12) 
at radii > 1". More flexible disk and/or bulge models are needed to flt this region 
satisfactorily. 

• We have not explored the effects of noise due to the limited number of bright stars 
in the nucleus. Even in the brightest parts, there are fewer than 10^ stars per (O'.'l)^ 
area (L98), suggesting that the uncertainties due to Poisson statistics are at least a few 
percent in the STIS kinematics; the errors are likely to be even larger in the photometry 
because deconvolution enhances Poisson noise. 

• We have fit the kinematics only along the P1-P2 axis. Integral- field spectroscopy 
by BOl provides the mean velocity and velocity dispersion across the entire two- 
dimensional face of the nucleus, although at lower spatial resolution. We have shown 
that this data is approximately consistent with our model (Figures 9 and 10) but it 
should be included in the model fits as well. 

A striking feature of our simulations is that the non-aligned model fits the data sub- 
stantially better than the ahgned model. As stated above, this conclusion is unhkely to be 
an artifact of our neglect of the self-gravity of the disk, since the improved fit is apparent 
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in both the photometry and the kinematics. Our best-fit inchnation is 54° compared to an 
inchnation of 77° for the large-scale M31 disk. This inclination is consistent with values 
derived by other investigators when modeling the nucleus as a thin disk: 55° (BOl), 50° 
(Peng 2002), and 52° (Sambhus & Sridhar 2002). The position angle of the line of nodes of 
the disk on the sky plane in the non-aligned model is 9i + 90° = 47°, which is reasonably 
close to the angle derived by BOl (54° for the model in their Fig. 20). 

Most of the fitted parameters in the aligned and non-aligned models are rather similar 
(Table 1). However, the non-aligned model has a lower peak in the mean eccentricity (Fig. 
1); this is turn is mostly controlled by the parameter a (cq. 12), which is 30% smaller in 
the non-aligned model. In the observations, a primarily affects the positions and amplitudes 
of the velocity maxima, while not greatly influencing the dispersion curve. Decreasing a 
in the aligned model significantly improves the fit to the velocity curve. The parameter oq, 
which controls the surface- brightness scale length of the disk (eq. 17), is more than a factor of 
three larger in the aligned model than in the non-aligned model. Decreasing in the aligned 
model brings more stars near the BH and therefore can dramatically increase the dispersion 
peak; this parameter does not greatly affect the rotation velocity. However, if these changes 
to a and are made to the aligned model, the fit to the photometry worsens dramatically, 
giving a P2 peak that is very much brighter than PI, and not matching the general features 
of the two-dimensional photometry. Conversely, if one changes the orientation angles of the 
non-aligned model to ahgn the disk with the large-scale disk of M31, both the photometric 
and kinematic fits become worse. 

Thus, it appears that an aligned eccentric disk can match the photometry or the kine- 
matics but not both simultaneously. 

There are some arguments that non-aligned nuclear disks may be present in other galax- 
ies, (i) The axis of the 0.2-pc radius mascr disk in the center of NGC 4258 is 119° from 
the axis of its host galaxy (Miyoshi et al. 1995). (ii) The SO galaxy NGC 3706 contains 
an edge-on disk of radius ~ 20 pc that is tipped by ~ 30° to the major axis of the ellip- 
tical isophotes at larger radii (Lauer et al. 2002). (iii) The jets in Seyfert galaxies, which 
are presumably perpendicular to the inner parts of the nuclear disk, are not perpendicular 
to the large-scale galaxy disk (Schmitt & Kinney 2002). A counter-argument is that the 
isophote twists that should be associated with non-aligned disks do not appear to be seen 
near the centers of edge-on disk galaxies, although this argument has not been quantified in 
a well-defined sample. 

The possibility that the nuclear disk is not aligned with the large-scale M31 disk raises 
several interesting issues. A natural question is whether the inner bulge of M31 shares the 
orientation of the nuclear disk or the large-scale disk. Normally, of course, the bulge of 
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M31 is assumed to be aligned with the large-scale disk, since bulges and disks appear to be 
aligned in edge-on spiral galaxies. A bulge aligned with the large-scale disk must be triaxial 
or barred (Stark 1977; Stark & Binney 1994; Berman & Laurent 2002). Ruiz (1976) has 
investigated the alternative possibility that the bulge is axisymmetric but not aligned with 
the large-scale disk; in this case she finds that fitting the photometric and kinematic data 
requires an inclination of ~ 55°. The close agreement of this inclination with the derived 
inclination of the nuclear disk is striking but inconclusive. 

If the nuclear disk is not aligned with the inner bulge, then it is subject to dynamical 
friction from the bulge, which usually damps its inclination relative to the bulge. Dynamical 
friction arises because of precession of the line of nodes of the nuclear disk on the equa- 
torial plane of the bulge. Numerical simulations of the damping of large-scale warps in 
non-spherical halos indicate that the damping time is typically not much longer than the 
precession time (Dubinski & Kuijken 1995; Nelson & Tremaine 1995). The precession time of 
the nuclear disk in the potential of the bulge is ~ 10^ yr, unless the bulge is very accurately 
spherical (the axis ratio of the isophotes of the inner bulge is 0.8-0.9; see Kent 1989 and 
Peng 2002). Thus we expect that a long-lived nuclear disk should be aligned with the bulge. 
A possible way to evade this conclusion would be to argue that dynamical interactions with 
the bulge excite, rather than damp, the inclination of the nuclear disk (Nelson & Tremaine 
1995); one attraction of this approach is that one might also find that interactions with the 
bulge excite the mean eccentricity of the disk. 

The current photometric and kinematic data on the nucleus of M31 are sufficiently 
accurate and feature- rich to justify the development of more accurate and fiexible dynamical 
models that incorporate the improvements listed at the beginning of this Section. Such 
models should yield new insights into the formation and structure of eccentric stellar disks. 
In the more distant future, we may look forward to high-resolution imaging of the nucleus 
using adaptive optics, interferometric imaging by the Space Interferometry Mission, and 
measurements of proper motions in the nucleus, by HST or its successors or by SIM. In 
this spirit. Figures 17 and 18 show the kinematic parameters of our models as viewed at a 
resolution of 0'.'02. 

We thank Ralf Bender, Karl Gebhardt, John Kormendy and Tod Lauer for expert 
advice, S. Sridhar for thoughtful discussions, and the referee, Eric Emscllem, for many 
suggestions that improved the presentation and content of the paper. We also thank Bender 
and his collaborators for providing their STIS observations of the M31 nucleus in advance of 
pubhcation. This research was supported in part by NSF grant AST-9900316 and by NASA 
grant HST-AR-09513.01-A. 
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Fig. 17. — The bulge-subtracted rotation speed (left) and velocity dispersion (right) as 
defined by the Gauss- Hermite expansion (19), as observed along a slit of width 0'.'02 at 
PA=39°. The non-aligned model is shown using a solid line and the aligned model using a 
dashed line. The curves are jagged because of Poisson noise in our Monte-Carlo simulation, 
which is smaller than the expected Poisson noise from the limited number of stars in the 
actual nucleus. 




Fig. 18. — The Gauss- Hermite parameters hs (left panel) and /i4 (right panel), as observed 
along a slit of width Qf!02 at PA=39°. The non-aligned model is shown using a sohd fine and 
the aligned model using a dashed line. 
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